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In this study we implemented a new imaging method to fuse functional near infrared 
spectroscopy (fNIRS) measurements and functional magnetic resonance imaging (fMRI) 
data to reveal the spatiotemporal dynamics of the hemodynamic responses with high 
spatiotemporal resolution across the brain. We evaluated this method using multimodal 
data acquired from human right finger tapping tasks. And we found the proposed method 
is able to clearly identify from the linked components of fMRI and fNIRS where and when 
the hemodynamic signals are changing. In particular, the estimated associations between 
fNIRS and fMRI will be displayed as time varying spatial fMRI maps along with the fNIRS 
time courses. In addition, the joint components between fMRI and fNIRS are combined 
together to generate full spatiotemporal "snapshots" and movies, which provides an 
excellent way to examine the dynamic interplay between hemodynamic fNIRS and fMRI 
measurements. 
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INTRODUCTION 

Functional near infrared spectroscopy (fNIRS) offers unsur- 
passed high temporal resolution and provides quantitative hemo- 
dynamic information for both oxyhemoglobin (HbC^) and 
deoxyhemoglobin (HbR), which plays an important role in the 
in vivo study of cognitive processing in the human brain (Jobsis, 
1977; Cope and Delpy, 1988; Hoshi, 2003; Singh et al., 2005; 
Huppert et al, 2009; Ye et al, 2009; Yuan et al, 2010a,b; 
Brunno et al., 2011; Egetemeir et al., 2011; Gagnon et al., 
2012; Yuan, in press). Similar to its fNIRS counterpart, func- 
tional magnetic resonance imaging (fMRI) is also a non- invasive 
imaging method that measures the hemodynamic responses 
to even-related neural activity with excellent spatial resolution 
and low temporal resolution. To take advantages of the com- 
plementary information from these two imaging modalities, 
two broad methods for fNIRS and fMRI/MRI integration have 
been developed for different clinical cases: (1) "spatial con- 
straint," in which spatial information from fMRI/MRI images 
are utilized to aid diffuse optical imaging based on fNIRS 
measurements (Carpenter et al., 2007; Ferradal et al., 2013); 
(2) "temporal correlation," where fMRI bold signals are pro- 
cessed to generate the correlation with Hb02 and HbR concen- 
tration changes converted from fNIRS recordings (Cui et al, 
2011; Tak et al, 2011; Gagnon et al, 2012). However, so far 
the linking between Hb02/HbR signals and fMRI spatial maps 
has not been extensively investigated, which represents one of 
the main challenges for fNIRS and fMRI fusion. Therefore, 
it is crucial for us to develop new imaging techniques to 
reveal the connections between these two measurements so that 
we are able to examine the dynamic interplay between space 
and time of hemodynamic responses with high spatiotemporal 
resolution. 



Interestingly joint independent component analysis (jICA) 
method has been developed to compute the linked temporally 
independent event-related potential (ERP) components and spa- 
tial independent fMRI components, which enables inferences to 
be made using estimated associations between fMRI sources and 
ERP electromagnetic sources (Calhoun et al., 2006a,b; Sui et al., 
2009). However, what has not been tried is a joint estimation of 
the temporal parts of fNIRS waveforms and the spatial maps of 
fMRI images. In this study, the jICA is extended to identify the 
spatiotemporal decompositions composed of fMRI spatial com- 
ponents indicating where the hemodynamic signals are changing 
and fNIRS components indicating when the hemodynamic sig- 
nals are changing. The fNlRS-fMRI fusion method will involve 
calculating for given stimuli, the connection between the time- 
locked fNIRS waveforms and fMRI activation maps for all par- 
ticipants or different measurement sections from a single subject. 
It is anticipated that the results derived from jICA could be visu- 
alized by computing spatiotemporal "snapshot," which provides 
an effective way to examine the dynamic interplay between fNIRS 
and fMRI hemodynamic sources. 

METHODS 

In terms of jICA, the joint spatial and temporal independences of 
fNRIS and fMRI are assumed to satisfy the following generative 
model for the data (Calhoun et al, 2006a,b), 

j^fNIRS _ fNIRS j^fMRI _ ^<jfMRI 

in which x^ 1 ^ is the group data from the chro- 
mophore concentration change of Hb02 or HbR for 
n subjectsM sections of a single subject and X^ 1118 = 

[xf< IR8 , IRS , ...X^ IRS f, X^ 1 is the group data 
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FIGURE 1 | (A) The data review for AHbC>2 measurements; (B) the data AHb02 results (middle) and AHbR results(right side). The axes (bottom 

review for AHbR measurements; (C) Channel configurations along the scalp middle and bottom right) illustrate the time scale, in ms, whereas the scale 
for right finger tapping tasks; (D) the section averaged fMRI image (left side), for the middle and right columns records chromophore concentrations in |xM. 
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FIGURE 2 | jICA decomposition of AHb0 2 and fMRI joint data for right 
finger tapping tasks: three components were found to be significantly 
correlated with finger tapping tasks [Independent components (ICs) 2, 
4, and 5] while two were well-correlated with the physiology and body 
movement noise (ICs 1 and 3). Each of these components is shown in a 
separate panel in the figure: IC2 (A), IC4 (B), IC5 (C), IC1 (D) and IC3 (E). 
The fMRI maps are thresholded at |Z| > 1 .5 for display purposes. 
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FIGURE 2 | Continued 

The averaged event-related AHbC>2 time course is shown in yellow (the 
same for all figures) and the AHb02 component is plotted in cyan. Positive 
(orange) and negative (blue) Z-values are shown in the image. The axes 
(bottom) illustrate the time scale, in ms, whereas the scale (left) records 
chromophore concentrations in |xM for the figures on the left column. The 
scale for the figures on the right column shows the bold signal intensity. 
LM, Left primary motor cortex; SMA, Supplementary motor area; R Parietal 
cortex; V, Visual cortex; BR Body movement/physiology noise. 
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It is note Equation 1 can be rewritten as a single matrix equation, 
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We employ the infomax ICA method for jICA of Equation 3, 
which utilizes a gradient ascent iteration algorithm to maximize 
the entropy of the output of a single layer neural network (Bell 
and Sejnowski, 1995). The resulting updated equation for the 
algorithm to calculate the shared unmixing matrix W (i.e., the 
inversion of A), the fused independent fNIRS sources u miRS and 
fMRI sources ifl^ 1 is as follows, 



Aw = )] u 



■ 2y mms (u mms f 



■ 2y mm (u MRl 



W (4) 



m wnicn y mRS = g(u mRS ), y 

WXfMRS^ M fMRi = wx MR \ and g(x) = 1/(1 + e~ x ) is the 
non-linearity in the neural network. The initial value for W, 
W(0) is a matrix composed of random vectors. 

The jICA method generally doesn't provide us the details on 
how the components between fNIRS and fMRI interact with one 
another. To achieve this, the spatiotemporal "snapshots" of the 
most significant components are generated in two ways (Calhoun 
et al, 2006a,b). First we need to calculate the linear combina- 
tion of the fMRI components weighted by their joint fNIRS time 
courses for a specific point in time. If the n spatial (fMRI) and 
temporal (fNIRS) components are given by S = [si 52 . . . s n ], 
and T = [t\ ... t n \ where t \ is a T x 1 vector containing 
the fNIRS time points and s; is a V x 1 vector represents the V 
brain voxels, the fMRI movie is computed as MfMRi = x S T . 
It is noted the absolute value is utilized here because the joint 
components are fused using a single parameter. And a change in 



„fMRI 



= g(u MR1 ), u mRS = 
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FIGURE 3 | jICA decomposition of AHbR and fMRI joint data for right 
finger tapping tasks: two components were found to be significantly 
correlated with right finger tapping task (ICs 1 and 3) while three were 
well-correlated with the physiology and body movement noise (ICs 2, 

4, and 5). Each of these components is shown in a separate panel in the 
figure: IC1 (A), IC3 (B), IC2 (C), IC4 (D) and IC5 (E). The fMRI maps are 
thresholded at |Z| > 1.5 for display purposes. The averaged event-related 
HbR time course is shown in yellow (same for all figures) and the AHbR 
component is plotted in cyan. Positive (orange) and negative (blue) Z-values 
are shown in the image. The axes (bottom) illustrate the time scale, in ms, 
whereas the scale (left) records chromophore concentrations in |xM for 

(Continued) 



FIGURE 3 | Continued 

the figures on the left column. The scale for the figures on the right 
column shows the bold signal intensity. LM, Left primary motor cortex; 
SMA, Supplementary motor area; R Parietal cortex; V, Visual cortex; BR 
Body movement/physiology noise. 



the amplitude of the fMRI component is directly linked to the 
change in the fNIRS component by this parameter. Meanwhile 
the fNIRS movie is estimated by Mf^iRs = T x |S| r , in which the 
time course for a given fMRI voxel is computed. 

RESULTS 

BEHAVIOR TASKS AND fNIRS-fMRI RECORDINGS 

The fNIRS tests are implemented with a block design for a right 
finger tapping task. The experiment is performed with a 24- 
chnnel fNIRS system (Oxymon MKIII, Artinis), which has 8 
sources, 8 detectors, and 24 channels. In this system, two con- 
tinuous wave lights at wavelengths 781 and 856 nm are emitted 
at each source fiber. In the case of block design for right fin- 
ger tapping tasks, the onset time for the first trigger was at 42 s, 
then followed by a 21s period of activation alternated with a 
30 s period of rest. This was repeated 10 times for the sub- 
ject. As such, the total recording time was 552 s. During the 
task period, subject was instructed to perform a finger flex- 
ion and extension action repeatedly. Data segmentation, which 
is also known as epoching in signal processing, is utilized to 
chop up the continuous fNIRS data into small time periods. 
The general way to do this is to extract segments surrounding 
the event codes from the experiments, e.g., from —35 s prior 
to the event onset until 35 s after the event code in this study. 
The original photon density datasets could be downloaded from 
(http://bisp.kaist.ac.kr/NIRS-SPM/Sample_data). The converted 
AHb02 and AHbR measurements (the unitless differential path 
length factor DPF = 4; sampling rate: 9.75 Hz) are filtered, seg- 
mented and plotted in Figures 1A,B, respectively. The configu- 
rations of 24 channels located on the scalp are also provided 
in Figure 1C. We found channels 8-12 are very unique because 
their locations are close to the left motor cortex so that they are 
able to reveal the hemodynamic responses from both the right 
finger tapping stimuli and other events. As a result, signal aver- 
aging is implemented for the datasets from these five channels 
to generate the representative AHb02l AHbR measurements with 
reduced noise. In particular, five-section (runs 2, 4, 6, 8, and 10) 
measurements are selected from the representative ten- run fNIRS 
recordings for further fusion analysis since they show the best 
signal noise ratio during the stimulus processing. 

For fMRI recordings, the fiber length was 10 m to connect the 
optodes in the MR scanner to the NIRS instrument in the MR 
control room. A 3.0 T MRI system (ISOL, Republic of Korea) was 
used to measure the BOLD response. During the blocked task 
paradigm, the echo planar imaging (EPI) sequence was used with 
TRITE = 3000/35 ms, flipangle = 80°, 35 slices, and 4 mm slice 
thickness (Ye et al., 2009). The fMRI data could be downloaded 
partly from the same website for the same subject or requested 
from the KAIST lab. Based on the datasets provided, we gener- 
ated five mean MRI images which are considered as five-section 
fMRI measurements (Calhoun et al, 2001; Ye et al, 2009). 
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FIGURE 4 I AHb02 components and fMRI "snapshots": on the left of 
each window is shown a linear combination of the fMRI maps that are 
well-correlated with finger tapping tasks (composite ICs 2, 4, and 5), 
weighted by the AHb0 2 part of the components at a specific point in 
time. On the right of each window is shown the estimated AHbC>2 
components that are correlated with right finger tapping tasks. The time 
courses for IC2 (in blue), IC4 (in green) and IC5 (in red) are also plotted on 
the right of each window. Such a display provides a dynamic way to 
visualize the brain activity at different time points: -20943 ms (A), 
-11217 ms (B), -8210 ms (C), -5559 ms (D), -2729 ms (E), 2044 ms (F) 
and 10884 ms (G). 



FIGURE 4 | Continued 



Then we will combine the five-section averaged fNIRS record- 
ings in temporal domain and five mean fMRI images in space 
to extract the joint independent components for hemodynamic 
fusion. It should be noted here the mean fMRI images are not 
generated from the raw fMRI images at each TR. As such, the 
fMRI image for each section may not accurately match its cor- 
related fNIRS measurements in terms of recording times, which 
may have some influence on the final results. However, the influ- 
ences should not be that significant since we use mean data from 
each section for fusion. The five-section averaged fMRI image and 
AHb02/ AHbR results are presented in Figure ID for references. 

RESULTS OF fNIRS-fMRI FUSION 

The computed joint components are provided in Figures 2, 3, 
in which the spatial components and regions of brain activity 
from the fMRI maps are plotted on the right while the correlated 
temporal joint components along with the mean AHb02l AHbR 
signals are given on the left. We can see from Figures 2, 3 that 
the time courses of temporal AHb02l 'AHbR components cor- 
respond well to different peaks on the mean AHbC>2l AHbR 
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FIGURE 5 | AHb02 components and fMRI "snapshots": on the left of 
each window is shown a linear combination of the fMRI maps that are 
well-correlated with physiology and body movement noise (composite 
ICs 1 and 3), weighted by the ts.Hb0 2 part of the components at a 
specific point in time. On the right of each window is shown the 
estimated AHbC>2 components that are correlated physiology and body 
movement noise. The time courses for IC1 (in blue) and IC3 (in green) are 
also plotted on the right of each window. Such a display provides a dynamic 
way to visualize how the noise affects the brain activity at different time 
points: -20584 ms (A), -12806 ms (B), 984 ms (C), 5757 ms (D) and 
11062 ms (E). 



signals. Interestingly, when put together the fMRI maps and 
fNIRS signals, we found the five spatial joint components cor- 
relate very well with the five temporal joint components. For 
example, the first negative peak of AHb02 in Figure 2A basi- 
cally identifies the physiology noise such as visual activity and 
eye blinks. Also visible for this waveform of the second tem- 
poral joint component in Figure 2A is a late positive peak. 
During the positive peak, the hemodynamic activity is mainly 
found in the left motor cortex. Then the joint components from 
AHb02/ AHbR and fMRI are combined together to generate full 
spatiotemporal movies to more clearly display the dynamic inter- 
play between the fNIRS and fMRI measurements. In particular, 
two movies are created to show the spatiotemporal dynamics 
of the fMRI and AHbC>2 fusion, in which Movie 1 displays the 
hemodynamic responses of right finger tapping events (compos- 
ite spatiotemporal joint components 2, 4, and 5) while Movie 2 
basically reveals brain activity from physiology and body move- 
ment(composite joint components 1 and 3). The spatiotemporal 
hemodynamic changes derived from fMRI and AHbR fusion 
are also recorded by Movies 3, 4 for right finger tapping tasks 
(composite joint component 1 and 3) and physiology and body 
movement noise (composite joint component 2, 4, and 5), 
respectively. 

Figure 4 shows seven "snapshots" snipped from Movie 1 pro- 
duced from fused AHbC>2 and fMRI data for finger tapping 
tasks. On the left of the "snapshots" a linear sum of the fMRI 
maps weighted by their respective AHb02 time courses are pro- 
vided while the mean AHb02 waveforms combined the estimated 
AHbC>2 components are plotted on the right of each "snap- 
shot." Figure 5 plots five "snapshots" captured from fused the 
fMRI and AHb02 data, which basically reveals the spatiotem- 
poral responses of the physiology and body movement noise. 
Likewise, the "snapshots" from Movie 3, 4 produced from fused 
fMRI and AHbR components are given in Figures 6, 7 for fin- 
ger tapping tasks and physiology and body movement noise, 
respectively. 

DISCUSSION 

We show, for the first time, a spatiotemporal reconstruction of 
the human brain's responses with high temporal and spatial res- 
olution by fusing together the fNIRS signals and fMRI images. 
The jICA is a very key and unique technique for connecting 
hemodynamic fNIRS temporal signal and hemodynamic fMRI 
spatial information according to the joint constraints of temporal 
and spatial independences. Such direct correlation further pro- 
vides a very valuable tool to localize non-invasively those high 
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FIGURE 6 | Continued 



FIGURE 6 | AHbR components and fMRI "snapshots": on the left of 
each window is shown a linear combination of the fMRI maps that are 
well-correlated with finger tapping tasks (composite ICs 1 and 3), 
weighted by the AHbR part of the components at a specific point in 
time. On the right of each window is shown the estimated AHbR 
components that are correlated with right finger tapping tasks. The time 
courses for IC1 (in blue) and IC3 (in green) are also plotted on the right of 
each window. Such a display provides a dynamic way to visualize the brain 
activity at different time points: -13460 ms (A), -10452 ms (B), -1082 ms 
(C), 1747 ms (D), 4045 ms (E), 9527 ms (F) and 13768 ms (G). 



resolution structures that underlines temporally well- resolved 
fNIRS responses. Meanwhile, deep brain structures, which are 
not readily detectable with scalp fNIRS alone, can now be iden- 
tified by the fusion technique when aided by the correlated spatial 
components from fMRI. 

It is observed from Figure 2 that the first and third AHb02~ 
fMRI joint components basically identifies body movement and 
negative eye blink activities in pre-frontal lobe while the second, 
fourth and fifth AHb02~fMRI joint components are fused with 
hemodynamic activity in left primary motor cortex (PMC), sup- 
plementary motor area (SMA), and motor association cortex such 
as parietal cortex. In particular, the fifth joint component basi- 
cally identifies the neural activity in the SMA and parietal cortex. 
These observations are further validated by Movies 1, 2, which 
show the dynamic interplay between space and time of AHbC>2- 
fMRI hemodynamic responses. It is noted Movie 1 displays the 
spatiotemporal changes of neural activation patterns for finger 
tapping tasks. Further, we can see from Figure 3 that the sec- 
ond, fourth, and fifth AHbR-fMRI joint components basically 
reveal body movement and other physiology noise while the first 
and third AHbR-fMRI joint components identify the patterns of 
brain activity in the left motor cortex and motor association cor- 
tex. In particular, the first AHbR-fMRI joint component reveals 
brain activity in SMA and parietal cortex. These spatiotemporal 
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FIGURE 7 | Continued 



FIGURE 7 | AHbR components and fMRI "snapshots": on the left of 
each window is shown a linear combination of the fMRI maps that are 
well-correlated with physiology and body movement noise (composite 
ICs 2, 4, and 5), weighted by the AHbR part of the components at a 
specific point in time. On the right of each window is shown the 
estimated AHbR components that are correlated physiology and body 
movement noise. The time courses for IC2 (in blue), IC4 (in green) and IC5 
(in red) are also plotted on the right of each window. Such a display provides 
a dynamic way to visualize how the noise affects the brain activity at 
different time points: -14518 ms (A), -9391 ms (B), -3734 ms (C), 
7052 ms (D) and 13946 ms (E). 



findings from fused AHbR-fMRI data are further processed to 
generate Movies 3, 4, in which we can see the hemodynamic 
changes for a finger tapping task with high temporal and spatial 
resolution, and we can also see clearly how physiology and body 
movement noise corresponds to the neural activity. Importantly, 
we can observe from Figures 2, 3 as well as Movies 1, 3 that strong 
brain activity mainly occurs in the left PMC, which validates 
stimuli in the right finger tapping task will yield the activations 
in the left PMC. 

Alternatively, we can see from Figure 4 the estimated AHb02 
time courses (i.e., linear combinations of the significant AHb02 
components) at specific positions in the cerebral cortex. The 
first "snapshot" in Figure 4A is associated with physiology 
and body movement noise while the second to fourth "snap- 
shots" show the AHb02 concentration change during onset 
and stimulus processing of right finger tapping, which are val- 
idated by the spatial maps of fMRI image on the left. We can 
further observe from Figure 4 that the fifth "snapshot" cap- 
tures the maximized spatiotemporal hemodynamic responses 
while the sixth and seventh ones show the decreased trends 
of neural activity. So if we examine the different "snapshots" 
in Figure 4, the spatiotemporal dynamics of the right fin- 
ger tapping are revealed. In addition, the negative and posi- 
tive peaks of the estimated AHbC>2 time courses in Figure 5 
directly correspond to the physiology and body movement noise 
though partly they capture the weak and negative hemody- 
namic responses of right finger tapping in the case of strong 
body movement. 

The AHbR-fMRI "snapshots" in Figure 6 identify SMA and 
parietal cortex with early peak hemodynamic responses. Later 
responses occurring in left PMC show peaks of the AHbR 
responses during stimulus processing. In particular, "snapshots" 
2-7 display the hemodynamic changes in detail with high spa- 
tiotemporal resolution for finger tapping task. Figure 7 provides 
the portions of estimated AHbR time courses at specific posi- 
tions of fMRI images, in which we observe the peaks basically 
corresponding to the activities of physiology noise, eye blinks, and 
body movement. 

In summary, we have used a combination of hemodynamic 
fNIRS and fMRI data to visualize the neural activation patterns 
with high spatiotemporal resolution. The fusion techniques have 
shown the potential by using the joint hemodynamic data to 
identify unique neural information that cannot be revealed in 
either technique alone. In particular, we implemented a joint 
decomposition of fNIRS and fMRI data, which was linked or 
fused by a common mixing parameter. The present method does 
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not involve the solution of inverse problem for diffuse optical 
imaging or involve the use of the threshold for the fMRI data. 
However, we do need to decompose the combined data into spe- 
cific components by using jICA, which are composed of fNIRS 
and fMRI portions. Separating the data into joint components 
will provide us a useful way to examine component specific dif- 
ferences among different patient groups or different stimulus 
tasks. 

We have to point out the present investigation has sig- 
nificant limitations though we are encouraged that the 
results we have found seem meaningful and valid. For 
example, we only validate data from a single subject even 
without accurate measurements from fMRI. In future we 
would like to examine our method to incorporate differ- 
ent subject groups and complex neural stimuli as well as 
multiple time points acquired from multiple channels of 
fNIRS systems. 
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Movie 1 | AHb02 and fMRI fusion movie for finger tapping tasks. 

Movie 2 | AHb02 and fMRI fusion movie for physiology and body 
movement noise. 
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